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Abstract.  A  special  approach  to  deal  with  elliptic  problems  with  singularities  is  intro¬ 
duced.  It  is  shown  that  this  approach,  to  be  called  the  auxiliary  mapping  technique, 
in  the  frame  of  the  p- version  of  the  finite  element  method  yields  an  exponential  rate  of 
convergence  at  virtually  no  extra  cost  and  from  the  view-point  of  implementation  it  is  the 
easist  and  the  cheapest  method  of  all.  It  is  also  shown  that  this  technique  can  be  used  for 
elliptic  problems  on  unbounded  domains  in  R 2  as  well. 
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1.  Introduction 


In  the  theory  and  practice  of  the  finite  elements  a  considerable  effort  has  been  made 
to  design  special  approaches  to  deal  with  problems  on  domains  with  comers  (especially, 
cracks)  and  with  problems  on  infinite  domains. 

In  the  case  of  nonsmooth  boundaries  the  following  approaches  are  the  most  typical: 

(IA)  Mesh  Refinement([4],[6],[7],[8],[9],[13],[14],[25]). 

(IB)  Use  of  Special  Elements([l],[2],[22],[26j). 

(I C)  Use  of  the  Enriched  (nonlocal)  Bases  Functions([16],[23]). 

In  the  case  of  the  infinite  domains,  the  widely  used  approaches  are: 

(2A)  Domain  Restriction  Procedure([15j,[17]). 

(21?)  The  Boundary  Element  Method  and  Combination  of  the  Finite  Element  Method  and 
the  Boundary  Element  Method([18],[21]). 

(2C)  Infinite  Element  Approach([28],[29j). 

The  method  (1C)  is  in  practice  almost  never  used  because  this  approach  destroys  the 
architecture  of  the  finite  element  method  program.  The  method  (IB)  is  to  devise  special 
finite  elements  in  which  the  approximation  mimics  the  singularity  in  elements  neighboring 
singular  points.  By  this  the  accuracy  is  increased,  nevertheless  the  rate  of  convergence  of 
the  h-version  is  not  improved.  By  far  the  most  effective  approach  for  handling  singularities 
is  the  method  (1A)  but  the  success  of  it  depends  on  the  proper  mesh  refinement. 

The  method  (2 C)  is  not  often  used  in  practice.  The  Boundary  Element  Method  handles 
well  the  problem  with  constant  coefficients  on  domains  with  bounded  boundary,  but  this 
method  has  difficulties  when  the  coefficients  of  the  equation  are  not  constant,  the  boundary 
is  unbounded,  and  is  ineffective  when  the  differential  equation  is  not  homogeneous.  The 
hybrid  method  between  the  Finite  Element  Method  and  the  Boundary  Element  Method 
is  sometimes  used  but  lies  in  the  h-version  and  it  has  the  same  character  and  difficulty 


as  (IS).  The  most  widely  used  method  in  practice  for  dealing  with  infinite  domain  is 
the  method  (2.4).  However  the  main  concern  in  this  approach  is  to  choose  the  trimeating 
domain  properly  to  get  a  desired  accuracy. 

There  are  three  versions  of  the  finite  element  method:  the  h-  version,  the  p- version,  »nd 
the  h-p  version.  The  h-ve rsion  is  the  standard  one,  where  the  degree  p  of  the  elements  is 
fixed,  usually  on  low  level,  typically  p  =  1,2,3  and  the  accuracy  is  achieved  by  properly 
refining  the  mesh.  The  p-version,  in  contrast,  fixes  the  mesh  and  achieved  the  accuracy 
by  increasing  the  degrees  p  of  the  elements  uniformly  or  selectively.  The  h-p  version  is 
the  combination  of  both.  In  this  paper,  we  are  mainly  concerned  with  the  p-version  of 
the  finite  element  method([8j,[10j).  We  will  see  that  the  p-version  in  many  cases  avoids 
naturally  some  difficulties  arising  in  the  h-ve  rsion  because  it  uses  elements  of  large  size. 
The  basic  aspects  of  the  finite  element  approach  are 

(а)  The  geometry,  where  the  domain  is  defined  and  partitioned  into  elements. 

(б)  Creation  of  elemental  stiffness  matrices  and  elemental  load  vectors. 

(c)  Assembling  elemental  stiffness  matrices  and  load  vectors  together  to  form  the  global 
stiffness  matrix  A  and  the  global  load  vector  b 

( d )  Solving  the  system  Ax  =  b 

(e)  Post  processing. 

The  elemental  stiffness  matrices  axe  based  on  the  (elemental)  mapping  of  the  element 
onto  the  standard  element  and  the  polynomial  shape  functions  on  the  standard  element. 
The  usual  (elemental)  mapping  is  typically  of  isoparametric  or  blending  type. 

In  this  paper,  addressing  two  dimensional  problems,  we  will  show  that  the  effective 
use  of  other  mapping  (for  example  conformal  mapping  or  quasi  conformal  mapping)  in 
conjunction  with  the  usual  (elemental)  mapping  in  the  frame  of  the  p-version  of  the  finite 
element  method  yields  an  exponential  rate  of  convergence  at  virtually  no  extra  cost.  This 
approach  can  be  directly  and  practically,  without  any  changes,  implemented  in  the  p- 
version  of  the  Finite  Element  Code  PROBE([24j). 
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The  outline  of  this  paper  is  as  follows:  In  section  2  the  generalized  elements  are  intro¬ 
duced  and  the  conditions  for  elemental  mappings  to  yield  conforming  elements  are  stated. 
Also  listed  are  some  preliminary  lemmas.  In  section  3  we  prove  that  the  p-version  of  the 
finite  element  method  leads  to  an  exponential  rate  of  convergence  when  the  exact  solution 
of  the  problem  is  analytic.  In  section  4  the  method  of  auxiliary  mappings  is  introduced.  It 
is  shown  that  one  can  obtain  an  exponential  rate  of  convergence  at  no  extra  cost  when  this 
technique  is  applied  to  Laplace  equations  on  domains  with  comers.  In  section  5  we  show 
that  the  mapping  technique  also  works  for  general  elliptic  boundary  value  problems.  In 
section  6  this  technique  is  applied  to  the  problems  on  the  exterior  of  the  bounded  domains. 
Finally,  the  conclusions  are  summarized  in  section  7. 

2.  Preliminaries 

Throughout  this  paper,  Q  denotes  a  simply  connected  (bounded  or  unbounded)  open 
subset  of  R 2  ,  where  R 2  denotes  the  usual  Euclidean  space  with  x  =  (ii ,  x2)  E  R2  .  As 
we  mentioned  above,  the  domain  Q  under  consideration  is  partitioned  into  the  elements 
e;,i  =  1,2,  ...n.  By  an  dement  e  we  mean  an  open  curvilinear  triangle  or  rectangle.  A 
curvilinear  triangle  is  a  domain  which  is  bounded  by  three  smooth( analytic)  arcs( sides) 
whose  ends  axe  called  the  vertices.  A  curvilinear  rectangle  has  analogous  meaning.  We 
allow  the  triangle  to  be  unbounded,  where  one  vertex  is  located  at  oo  analogously  as  in 
the  theory  of  complex  variable.  Then  S7  =  U,e7  ,  where  by  the  bar  we  mean  the  closure  in 
R2  .  As  usual,  we  will  assume  that  e<  D  e;  is  either  empty  or  is  a  common  vertex  or  is  an 
entire  common  side. 

By  T  and  Q  we  denote  the  standard  triangular  element  and  the  standard  quadrilateral 
element  respectively  as  shown  in  Fig. 2. 1(a)(6).  The  coordinates  in  the  standard  plane  are 
denoted  by  £  =  (fi,^)-  As  usual  we  define  on  the  element  e  the  finite  element  space  5(e) 


which  is  the  span  of  elemental  shape  functions,  denoted  by  Nf(x),x  6  e,  i.e.  5(e)  = 
span,yVte(x). 

By  the  standard  shape  function  Wi(f)  we  mean  the  shape  functions  on  the  standard 
element  E(E  =  T  or  E  =  Q)  which  are  divided  into  the  groups  of  (1)  nodal,  (2)  side 
and  (3)  internal  shape  functions.  We  assume  as  usual  that  a  bijective  smooth  elemental 
mapping  from  E  to  e  is  such  that  $e(£)  =  x  €  e,  for  f  €  E,  where  E  =  T  or  E  =  Q 
according  as  e  is  (curvilinear)  triangular  element  or  (curvilinear)  quadrilateral  element. 

Then  the  elemental  shape  functions  Nf  are  related  to  the  standard  shape  function  by 
Nf(x)  =  Wi($71(;r))(  equivalently  by  N*($e(£))  =  Wi( £)).  Alternatively  we  denote  it  by 
N*  =  Ar,  o  $“*(  or  =  iV,e  o  $e).  We  associated  to  every  element  e,  of  the  partition  of  Cl 
the  mapping  •  We  will  assume  that  these  mappings  satisfy  the  usual  conditions  used  in 
the  finite  element  method  that  lead  to  conforming  finite  elements.  This  means  that  if  A  is 
a  common  vertex  of  the  element  ,  j  =  1,2,  ...m  and  N^’  are  the  nodal  shape  functions 
associated  to  the  vertex  A  then  the  function  v  such  that 

f  u(i)  =  N-."’  (x)  on  ek.,j  =  1,2,  ...,m 
(  v(x)  =0  for  x  £  eicj,j  =  1,2,  ...m 

is  a  continuous  function  on  Q  .  Similarly  if  7  is  the  common  side  of  elements  etj  and  ejt2 
and 

=  1,2 

are  the  side  shape  functions  associated  to  the  common  side  then  the  function  v  defined  by 

(v(x)  =W*‘>( x)  on  ekj,j  =  1,2. 

\v(x)  =0  for  x  $  eknj  =  1,2. 

is  a  continuous  function  on  Cl  .  Finally  the  supports  of  the  internal  shape  functions  are  in 
the  element. 
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There  axe  many  possible  elemental  mappings.  We  will  impose  some  restriction  on  the 
elemental  mappings.  Especially  we  will  assume  that  $e  has  the  circular  arc  property:  Let 
e  be  a  curvilinear  triangle  with  a  circular  side  as  shown  in  Fig.  2.2,  then  $e  linearly 
transforms  the  length  of  arc  of  the  standard  triangle  on  the  length  of  arc  of  the  element 
e.  If  e  is  a  curvilinear  quadrilateral  element  with  a  circular  side,  the  meaning  of  3>e 
with  circular  property  is  analogous.  Let  us  note  that  the  code  PROBE([24])  is  using  the 
elemental  mappings  with  circular  arc  property. 

In  general  we  will  consider  two  domains  ft  and  Cl  in  R2  such  that  'P(fl)  =  Cl  where  'P  is 
a  smooth  bijective  mapping  of  ft  onto  ft  .  If  u  is  a  funtion  defined  on  ft  then  u  will  denote 
the  function  defined  on  ft  by 

u  =  uo'P~1(oru  =  uo'P) 

Obviously  ifft  =  e,ft  =  T(orQ)  then  $  =  as  a  special  case  of  notation  we  have  used 
above.  When  there  is  a  mapping  <p  which  transforms  a  curvilinear  triangular  (quadrilateral) 
element  e  onto  an  element  e*  of  the  same  type,  we  will  construct  the  elemental  mapping 
<Pe-  of  the  standard  element  E(E  =  T  or  Q)  onto  e*  so  that 
M  o  <p-‘  =Nf  =(iV,e)  09-1 

=  JV,  0  o  for  any  shape  function  W,  on  E. 

Note  that  if  <p-1  is  a  mapping,  from  an  element  e"  with  a  circular  side  onto  another  element 
e  with  a  circular  side,  which  is  linear  in  the  arc-length,  and  if  <Pe.  has  circular  arc  property 
then  the  composite  mapping  <^_1  0  <Pe.  also  has  the  circular  arc  property.  We  have  seen 
that  the  main  objective  is  to  create  the  set  5(e)  of  elemental  shape  functions  which  led 
to  the  conforming  base  functions.  We  can  now  use  not  only  one  elemental  mapping  but  a 
family  of  such  mappings.  For  example,  let  >p>~1  maps  e*,  an  element  with  circular  side, onto 
e  .  an  element  with  circular  side,  then  we  can  enrich  the  space  of  elemental  shape  functions 
for  e  by  functions  W,  0  ^r.1  0  <jj,  where  W,  denotes  the  standard  shape  functions  which  are 
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Fig.  2.1.  (a)  The  Standard  Triangular  Element,  (b)  The  Standard 

Quadrilatral  Element. 


% 

©  - 

Fig.  2.2.  An  Elemental  Mapping  of  the  Standard  Triangular 
Element  T  onto  a  Curvilinear  Tringular  Element  e. 
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zero  on  the  side  of  the  standard  element  that  is  not  associated  to  the  circular  side.  So 
far  we  mentioned  only  the  objective  to  construct  the  elemental  shape  functions.  The  next 
objective  is  to  construct  them  so  that  local  stiffness  matrices  can  be  easily  computed. 

Returning  now  to  the  general  mapping  'I'  of  ft  onto  Cl  ,  we  denote  the  Jacobian  of  'P~1 
by  Jr('I'-1)  and  its  determinant  by  |  Jr('I,~1)|  and  we  will  assume  it  is  positive  on  Cl  .  Then 
we  define 

M  =  (|J($-1)|)1/2(J('I')otf-1) 

and  we  get  an  obvious  lemma. 


Lemma  2.1.  Suppose  A  =  [a,j],  ai}  6  L^Cl),  1  <  i,  j  <  2  and  if  Q  =  M {ai;  oT  1  }Mt  = 
{  qtJ } ,  then  we  get 


(2.1) 


HXa'i(x) 


du  dv 
dx,  dxj 


du  dv 

WkWi 


and  quit)  €  L  OO  (Cl). 


We  are  particularly  interested  in  the  special  case  when  A  =  I, the  identity  matrix  and 
is  a  conformal  mapping  of  Cl  onto  Cl  .  Then  we  have 


Corollary  2.1.  If  'i  is  a  conformal  mapping  and  A  =  /  then  Q  =  I,  that  is, 


JJvu.  Vr dx  =  JJ  Vu  •  Vvd£. 


Corollary  2.2.  Let  A  =  {atJ(i)}  be  a  symmetric  2x2  matrix  such  that 
|aj>(x)|  <  7i  <  oo,  «t}(x)Thri}  >  7  2(i)]  +  qj),  72  >  0, 
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for  all  ( rjj ,  rj2 )  6  R2  and  if  $  is  a  conformal  mapping,  then  we  have 


qkliO  <  7*  <  X>i(Owi  -  72^i  +  *72) .72  >  Oi 

for  any  (rp,  rj2)  €  R2  ,  where  7*  and  72  depend  on  the  mapping  but  not  on  £  .  Hence  the 
conformal  mapping  preserves  the  eilipticity  of  the  operator 


Let  us  note  that  Corollary  2.2  is  also  true  under  a  weaker  condition  on  'F  ,  for  example, 
when 

=  |J(^)|diag(d1(x),d2(x)) 

and 

0  <  ki  <  di(x)  <  k2,  for  all  x  6  fh 

Finally  we  would  like  to  mention  a  conformal  mapping  which  will  play  an  important  role 
in  forthcomming  sections: 

z  =  £Q ,  2  =  X!  +  z'x2, (xj , x2)  (z  e(  resp  fl) 

S  =  6  +  *6,(6, £2)  €  e*(  resp  fi) 

and  this  mapping  characterized  by  the  constant  a  will  be  denoted  by  ipa  (  resp  $“)  re¬ 
spectively. 
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3.  The  p- Version  of  the  Finite  Element  Method 

Let  fl  C  R2  be  a  bounded  domain  whose  boundary  T  is  composed  of  smooth  (analytic) 
arcs  T,,  i.e.  T  =  U,r,. 

Consider  now  the  following  model  problem 


(3.1a) 


— Au  =  /  in  fl, 


(3.1b) 


u  =  0  on  Tq  =  Ui62>r,, 


(3.1c) 


onF/v  =  Ut6^rt, 


where  V  fl  Af  =  d>,  r.v  U  To  =  T,  and  g  £  ^(Fn).  For  simplicity  we  assume  that  T  d  ^ 
(otherwise  the  solution  would  not  be  unique  and  we  would  have  to  use  the  standard 
approach  which  would  avoid  this  technical  difficulty). 

Let  Hx(Cl)  be  the  usual  Sobolev  space  and  HlD  =  {u  £  Hl(Q.)  :  u  =  0  on  I'd}  then  as 
usual  we  will  consider  the  weak  solution  of  the  problem:  Find  u0  €  such  that 


where 


B(u0,v)  =  F(v),  for  all  v  £  Hq(CI) 


Vu  ■  Vvdx, 


B(u'v]= SL 

F(v)  =  JJ  fvdx  +  J 


fvdx  +  /  gvds. 
r* 


Furthermore  we  denote  ||u|| ^  =  B(u,u)  and  will  call  ||.(|e  energy  norm.  Because  of  the 
assumption  on  Fp  we  have  ||u||£  =  j| tz || ( q ) . 

Let  us  now  consider  a  (fixed)  partitioning  of  fl  into  the  elements  e,,  i.e.  fl  =  Uef  as 
in  the  previous  sections.  Because  of  the  assumption  about  e,  the  vertices  of  fl  are  the 
vertices  of  some  elements.  We  will  assume  that  the  partition  and  the  (elemental)  mapping 
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$e  introduced  in  the  previous  section  satisfy  the  (technical)  conditions  listed  in  ([4],  [14]) 
that  are  satisfied  in  the  code  PROBE. 

Let 

Sp  =  {u  6  Hq(SI)  :  ufe  o  is  a  polynomial  of  degree  p  on  E  for  all  element  e}, 

where  E  is  T  or  Q  according  as  e  is  a  triangular  element  or  a  rectangular  element,  and  let 
Np  be  the  dimension  of  Sp.  Then  the  p- version  of  the  finite  element  method  is  now  to  find 
up  €  Sp  such  that 

B(up,  v )  =  F(v),  for  all  v  €  Sp 

and  then 

(3.2)  llup  ~  uojjs  =  min  j|u>  -  u0\Ie- 
We  have  now 

Theorem  3.1.  Let  u0  €  #d(^)  be  such  that 

(3.3)  lluo||//*(n)  S  CDkk\,k  =  1,2,... 

where  Hk(Ll)  is  the  usual  Sobolev  space  and  C  and  D  are  constants  independent  of  k,  then 

(3.4)  ||u0  -  a, ||e  <  Ce-^’ 
where  c  and  7  >  0  are  independent  of  p. 

Proof.  By  modifying  the  proof  given  in  [4],  one  can  easily  proceed  the  proof  of  this 
theorem.  Hence  we  will  only  outline  the  key  steps  of  the  proof. 

(i)  Denoting  u«  =  u|e  the  restriction  of  u  on  the  element  e,  we  define  ue  =  ut  o$e,  where 
is  the  elemental  mapping  of  the  standard  element  E  onto  e. 
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(ii)  We  construct  ut^p  £  H1(E)  which  is  a  polynomial  of  degree  p  on  E  and  approximates 
the  function  ue  in  H1(E).  Using  (3.3)  we  get 

(3-5)  \\ue-ueJHHE)<Ce-^P 

(iii)  Let  utp  —  ue<p  o  $“1,  then  ue<p  approximates  ue  on  e  and  satisfies 

l|ue  ~  ut'P\\HHe)  <  Cie_72P 

but  ue  p  <£ 

(iv)  We  construct  a  correction  term  on  every  e  so  that  the  correction  function  up  can  be 
-  member  of  HxD(fl)  and  satisfy  (3.4).  (QED) 

We  also  have 

Theorem  IS. 2.  Let  ft  be  a  polygonal  domain  and  ifuo  £  Hq(fl),  q  >  1  integer  or  fractional, 
then  we  have 

If  uq  =  crP<p(9)x(r,  9),  where  ( r,d )  are  the  polar  coordinates  at  the  vertex  A  of  Q  ,<p(9)  is 
a  smooth  function  and  x(ri@)  JS  3  smooth  cut-off  function,  then 

ll»o  -  U,lle  <  CN-0 

For  the  proof,  see  [8].  (QED) 
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Example  3.1  Let  Q  be  an  L-shaped  domain  partitioned  into  12  elements  as  shown  in 
Fig.  3.1  and 

'D  =  {1,2} 

N  =  {3, 4, 5, 6, 7, 8} 
g  =  1  on  I\,  i  =  4, 5,6, 7 
g  =  0  on  r i ,  i  =  3, 8 

Let  us  assume  /  =  0  and  denotes  the  circular  subset  of  0,  which  is  composed  of  the 
curved  triangular  elements  in  Fig  3.1.  Then  on  Q  =  Q  —  the  solution  u0  satisfies 
(3.3),  but  on  Q  it  only  satisfies  uq  €  fl),  where  e  is  an  arbitrary  positive  number. 

Moreover  the  solution  has  the  leading  term  r2/,3sin(2/3)0  near  the  re-entrant  comer. 
Hence  applying  theorem  3.2  we  get 

Up)  =  ll-ii~Y^~  £  Cl v-2'3 

\Wo\\e 

Fig.  4.2  shows  £o  in  dependence  on  Np  in  log  x  log  scale.  It  also  shows  the  slope  of  the 
rate  Np2^3 .  If  we  would  solve  the  problem  on  Q  then  we  would  get 

||«o-u,||*<Ce-V^ 

that  is,  the  exponential  rate.  Hence  we  see  that  the  error  decrease  algebraically  only  be¬ 
cause  the  function  from  Sp  are  not  able  to  approximate  well  the  solution  in  f Ir0. 
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4.  The  Method  of  Auxiliary  Mapping  for  the  Laplace  Operator. 

So  far  we  constructed  the  spaces  5(e)  by  using  the  standard  elemental  mappings  <$e 
which  were  based  on  the  blending  mapping  technique. 

Let  us  now  propose  the  question  about  optimal  mapping  technique  which  is  adjusted 
according  to  the  character  of  the  solution.  To  each  element  e  we  associate  a  mapping 
(see  Section  2)  which  maps  e  onto  e*  =  <^>c(e),  and  then  use  the  composite  mapping 

=  (<^e)_1  o 

for  constructing  5(e).  We  will  say  that  the  set  of  mappings  <pe  is  admissible  if  the  mappings 
{(ipe)-1  o  <J>e. }  will  lead  to  conforming  elements  (see  Section  2).  Now  the  problem  is 
following:  Describe  the  families  of  admissible  mappings  and  select  one  family  that  leads 
to  optimal  (in  some  sense)  results. 

Let  us  show  such  a  family  on  the  mesh  shown  in  example  3.1.  Let  us  map  the  element 
ej,j  =  7,  ..,12  (see  Fig.  4.1)  onto  the  element  e*  by  the  conformal  map  (<pa)~l,  where 
Kpa  :  =  z.  The  mappings  {(<p“)-1  '  j  —  7, ...,  12}  stemming  from  the  conformal  mapping 

(<pa)-1,a  >  3/2  and  the  identity  mappings  <pi,  i  =  1, 2,  ..6,  lead  to  the  conforming  elements. 
In  fact  (because  Ro  =  1)  the  curvilinear  triangles  e;,j  =  7,...,  12  are  mapped  by  the 
mapping  (<p")_1  onto  the  curvilinear  triangle  e*  and  the  mapping  is  linear  on  the  circular 
arc.  This  together  with  the  fact  that  4>e.  has  the  circular  arc  property  shows  that  elements 
constructed  by  <p"  o  are  conforming.  Using  Corollary  2.1,  we  see  that  in  our  example 
the  elemental  stiffness  matrices  of  ej,j  =  7, ...,  12  obtained  by  using  the  elemental  mapping 
$e>  =  (£>“  o  $e*  is  identical  with  the  elemental  stiffness  matrix  of  the  element  e*  obtained 
by  using  the  element  mapping  $e».  These  stiffness  matrices  are  directly  computable,  for 
example,  by  the  code  PROBE.  The  conformal  mapping  then  has  the  obvious  form: 
let  (r,  6)  be  the  polar  coordinates  in  the  z-plane  and  (p,  t/>)  the  polar  coordinates  in  the 

plane  ,  then  we  have 

r  =  pa  and  0  —  rpa\  p  =  r1/,Q!  and  tp  =  9/a 
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and  if  e  =  {(r,0)|O  <  r  <  Rq,(3q  <  6  <  (3\ }  then  e*  =  {(p,j/>)|0  <  r  <  do/a  <  r/>  < 

dj/a}.  Hence  we  compute  the  elemental  stiffness  matrix  for  e  in  e*  instead  of  in  e,  that 
is,  in  the  <f -plane. 

Assume  now  that  /  =  0  in  (3.1)  then  the  exact  solution  uo  can  be  written  in  the  form 
of  a  series 

oo 

(4.1)  u0(r,  9)  =  ^2  a*r(2/3)i:sin(2/3)&0  onQfio 

k=  1 

The  series  converges  absolutely  on  Qr0.  Now  using  the  conformal  mapping  <p°  :  z  = 
we  have  u0  =  Uo  o  <pa  and 

OO 

(4.2)  u0(p,i >)  =  Y^akp{2/3)ak  sin(2/3)t/>fca  on 

*= l 

Hence  if  a  =  (3/2 )m,m  an  integer,  uo  is  an  analytic  function.  Using  now  this  mapping  in 
the  finite  element  method  we  approximate  the  function  uo  on  e*  instead  of  the  function 
u0  on  tj.  Denoting  the  finite  element  solution  which  is  based  on  the  elemental  mapping 
pa  o  $e*  by  u",  we  get 

Theorem  4.1.  Let  us  consider  the  problem  given  in  example  3.1  with  f  =  0.  Suppose 
a  =  (3/2)m  ,  rn  is  an  integer.  Then  we  get 

K-Uo| \E<Ce~^’ 

where  C  and  7  are  independent  of  p  (and  hence  Np) 

Proof. 

Ilup  “  u0  ||  H1  (e  j  )  <  C\\u*  -  U0||//»(«;) 

and  apply  theorem  3.1  with  minor  modification.  (QED) 
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Tables  4.1  shows  the  relative  errors  for  m  =  2/3,  1,2  and  4.  Note  that  m  =  2/3  means 
the  case  when  no  mapping  technique  was  used.  Figs.  4.2  and  4.3  show  the  error  s  in 
log  5  x  log  .V  scale  and  Fig.  4.4(a)  shows  the  error  in  logs  x  %/N  scale. 

We  see  the  large  improvement  by  this  approach.  The  values  a  —  m( 3/2)  is  obviously 
the  most  advantageous  in  this  case.  Nevertheless  for  general  a>  l,a/  m( 3/2)  we  get  an 
asymptotic  improvement.  By  applying  Theorem  3.2,  we  get 


Figure  4.1. The  Scheme  of  and 


TABLE  4.1.  The  Relative  Error  for  a  =  ro(-),  m  =  2/3, 1,  2,  4 

when  /  =  0. 


DOF 


0.1176835578896779D  +  00 
0.2709809653Q06089D  -  01 
0.1392839460138052D  -  01  60 

0.4078061954970974D  -  02  100 

0.1137579818794599D  -  02  152 

0.3390584498845191D  -  03  216 

0.9705280386932562D  -  04  292 

0. 2757048098466325 D  -  04  380 


DOF 


1 

0.1628469724675572D 

+  00 

0.6774829063879436D 

-01 

0.4267764860979481D 

-01 

0.2930203635221640D 

-  01 

0.2236678198833828D 

-01 

0.17942319331 16408D 

-01 

0.1486700458517911D 

-01 

0.1261468544746016D 

-01 

0.2437842370603280D  +  01 
0.2757060901917393D  -  01 
0.1408437528088371D  -  01 
0.4142402706997820D  -  02 
0.1 149455073984022D  -  02 
0.3425494364421739D  -  03 
0.9736221988286052D  -  04 
0.2746260993089290D  -  04 


0.3913546088255135D  + 00 
0.1304975484223534D  +  00 
0.2328733675022739D  -  01 
0.4295634650207781D  -  02 
0.1237407887180692D  -  02 
0.4100441 132194179D  —  03 
0.1428178128929969D  —  03 
0.4942584844978042D  -  04 


Theorem  4.2. 

(4-4)  IK-u0|U< 

JVp 

where  C  is  independent  of  Np  . 

In  theorem  4.1  and  4.2  the  constant  C  depends  on  a.  We  have  seen  it  in  Table  4.1  and 
Fig.  4.2,  where  the  error  when  m  =  4  was  higher  than  that  when  m  =  1. 

To  see  more  in  details  these  effects,  let  us  investigate  the  error  of  an  approximation  of 
(x  +  l)a, x  €  I  =  (  —  1, 1)  by  polynomials  in  norm. 
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From  part  1  of  [12],  we  have 


Theorem  4.3. 

(4.5)  £p  =  Co(ot)  ^^(l  +  0{l/p))  as  p  — »  oo 
where 

co  (a) 

The  value  a  =  1/6  gives  the  rate  ep  =  0(p~4/3)  which  is  the  same  (in  p)  as  in  our 
example  when  no  mapping  is  used.  In  the  table  4.2  we  give  values  of  ep  for  some  values  of 
a  and  p  when  the  term  (1  +  0(l/p))  is  neglected.  We  also  give  the  values  of  the  norm  of 
(x  +  1)“  in  £2(-0  which  is  of  course  an  upper  bound  of  the  error.  We  selected  the  value  of 
a  so  that  |sina7r|  has  same  value  in  the  formula  (4.5).  We  can  see  the  typical  behavior, 
that  we  have  observed  in  Table  4.1,  namely  that  for  higher  a  the  results  are  worse  for 
low  p  (and  low  accuracy)  while  they  are  better  for  higher  degree  of  polynomial! and  higher 
accuracy).  Hence  in  general  the  choice  of  a  depends  on  the  required  accuracy. 

In  the  above  analysis,  we  have  considered  the  mesh  shown  in  Fig.  3.1.  We  can  of  course 
use  the  meshes  with  concentric  layers  which  are  used  in  the  h-p  version  too  (see  Fig.  4.5) 
and  use  the  mapping  as  before.  By  this  we  can  combine  mesh  refinement  approach  and 
the  auxiliary  mapping  technique. 

So  far  we  have  assumed  that  /  =  0  and  hence  the  solution  in  a  neighborhood  of  the 
corner  can  be  expressed  by  the  series  (4.1).  On  the  other  hand,  if  /  0,  /  G  Hs(ri)  then 

(4.1)  is  no  longer  true  and  the  solution  uq  in  f Irq  can  now  be  written  in  the  form(see,[19]) 

n  ^  n/  3 

(4.6)  u0(r,  0)  =  ^  a/tr(  sin(2/3)fcfl  +  ^  logrsin2M  +  v{r,d), 

k=o  k= 1 


T(1  +  a)2 1  sin7ra| 
(V^TT)7T 


where  v  G  H3+2(Qr0),  provided  that  (2/3 )(n  +  1)  >  s  +  1. 


Let  us  consider  /  =  1,  then  first  term  in  the  expansion  of  v  is  r2  log  r  sin  2kB  and  hence 
in  contrast  to  the  previous  case  we  can  not  achieve  (in  general)  an  exponential  rate  of 
convergence  even  when  /  is  analytic,  but  we  can  only  have  the  rate  N~2a  log  JVp  ([8], [9]). 
In  general  we  expect  that  the  positive  effect  by  the  auxiliary  mapping  technique  will  be 
slightly  smaller  for  /  ^  0  than  for  /  =  0.  Table  4.3  shows  the  errors  for  the  mapping 
when  a  =  m(§),  m  =  2/3, 1,2,  and  4,  and  m  =  2/3  stands  for  the  case  when  no  mapping 
technique  was  used  .  Fig.  4.6,  4.4(b)  show  the  behavior  of  the  error  in  the  logs  x  logiV 
scale  and  log  s  x  \/N  respectively.  We  actually  see  the  differences  as  predicted.  We  see 
that  the  auxiliary  mapping  technique  improves  the  accuracy  very  much  when  the  degree 
of  polynomial  is  >  4. 

So  far  we  have  used  only  one  mapping.  We  can  of  course  enrich  the  space  S(e)  of 
elemental  shape  functions  by  functions  stemming  from  identity  mapping,  i.e.  by  the  usual 
shape  functions.  We  can  also  add  analogous  shape  functions  using  the  conformal  mapping 
£  =  log  z.  By  such  an  enrichment  we  can  once  more  achieve  that  the  rate  of  convergence 
is  exponential  provided  that  /  is  analytic  in  Cl.  Of  course  the  complexity  in  computation 
now  is  higher  and  practical  questions  arise.  However  in  contrast  to  the  classical  enrichment 
approach  one  advantage  of  this  is  that  it  does  not  violate  the  finite  element  architecture(in 
the  implementation  we  add  only  internal  shape  functions  of  special  type  constructed  by 
the  mapping  procedure). 

As  we  mentioned  above  selecting  the  partial  ir  a  properly  and  using  higher  degree 
of  polynomials  yields  high  accuracy.  Let  us  also  state  that  we  can  use  high  degree  of 
polynomial  only  in  the  elements  with  vertex  at  the  corners  (in  our  example,  elements  7-12) 
while  keeping  a  low  degree  of  polynomials  in  ah  other  elements.  A  combination  with  mesh 
refinement  procedure  is  also  a  valuable  alternative  especially  because  of  the  easiness  of 
computing  the  local  stiffness  matrices. 
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TABLE  4.2.  The  error  ev  and  the  L2-norm  ||u||  of  u  =  (1  -f  x)a  on  I 

for  a  =  (n  +  1/6)  and  (n  +  5/6),  n  =  0, 1, 2,  3, 4. 


1 

7.4731D  -  2 

3.4188D  -  2 

3.2166D  -  2 

4.3432D  —  2 

5.9688D  -  2 

2 

4.3522D  -  2 

1.1596D  -2 

8.3258D  -  3 

6.5471D  -  3 

6.8665D  -  3 

3 

2.9657D  -  2 

5.3843D  -  3 

3.1913D  -  3 

1.7100D  -  3 

1.4805D  -  3 

4 

2.2025D  -  2 

2.9696D  -  3 

1.5168D  -  3 

6.0361D  -  4 

4.5034D  -  4 

5 

1.7272D  -  2 

1.8262D  -  3 

8.2602D  -  4 

2.5778D  -  4 

1.7031D  —  4 

6 

1.4063D  -  2 

1.2107D  -  3 

4.9412D  —  4 

1.2555D-4 

7.4850D  -  5 

7 

1.1769D  -  2 

8.4798D  -  4 

3.1661D  -  4 

6.7328D  -  5 

3.6720D  -  5 

8 

1.0059D  -  2 

6.1941D  —  4 

2.1381D  -  4 

3.8859D  -  5 

1.9592D  —  5 

9 

8.7406D  -  3 

4.6769D  -  4 

1.5049D  —  4 

2.3766D  -  5 

1.1170D  —  5 

10 

7.6975D  -  3 

3.6272D  -  4 

1.0953D  -  4 

1.5233D  -  5 

6.7187D  -  6 

11 

6.8543D  -  3 

2.8761D  —  4 

8.1951D  -  5 

1.0149D  -  5 

4.2242D  -  6 

12 

6.1605D  —  3 

2.3233D  -  4 

6.2760D  -  5 

6.9858D  -  6 

2.7565D  -  6 

13 

5.5809D  -  3 

1.9067D  -  4 

4.9023D  -  5 

4.9433D  -  6 

1.8565D  -  6 

14 

5.0904D  -  3 

1.5863D  -  4 

3.8951D  -  5 

3.5825D  -  6 

1.2850D  -  6 

15 

4.6707D  -  3 

1.3355D  -  4 

3.1412D  —  5 

2.6509D  -  6 

9.1077D  —  7 

IMI 

1.3747D  +  0 

1.5431D  +  0 

3.0238D  +  0 

2.3329D  +  0 

2.7495D  +  0 

1 

1.4586D  -  1 

2.5522D  -  1 

9.3989D  -  1 

1.9638D  +  0 

9.8958D  +  0 

2 

9.7720D  -  3 

1.3049D  -  2 

2.7987D  -  2 

4.4624D  -  2 

1.3096D  -  1 

3 

1.4357D  -  3 

1.5826D  -  3 

2.3129D  -  3 

3.0442D  -  3 

6.0879D  -  3 

4 

3.2433D  -  4 

3.0810D  -  4 

3.3440D  -  4 

3.7930D  -  4 

5.6332D  -  4 

5 

9.6187D  —  5 

8.0914D  -  5 

6.8870D  -  5 

6.9177D  -  5 

8.0567D  -  5 

6 

3.4420D  -  5 

2.6126D  -  5 

1.8105D  -5 

1.6411D-5 

1.5562D  -  5 

7 

1.4132D  —  5 

9.8131D  -  6 

5.6915D  —  6 

4.7192D  -  6 

3.7452D  -  6 

8 

6.4443D  -  6 

4.1370D  -  6 

2.0507D  -  6 

1.5720D  -  6 

1.0662D  —  6 

mm 

3.1925D  —  6 

1.9104D  -  6 

8.2288D  -  7 

5.8800D  -  7 

3.4655D  -  7 

1.6911D  —  6 

9.4970D  -  7 

3.6025D  -  7 

2.4157D  -  7 

1.2539D  -  7 

wM 

9.4678D  -  7 

5.0173D  -  7 

1.6947D  -  7 

1.0724D  —  7 

4.9565D  -  8 

12 

5.5527D  -  7 

2.7896D  -  7 

8.4689D  —  8 

5.0804D  -  8 

2.1105D  -  8 

13 

3.3880D  -  7 

1.6200D  -  7 

4.4555D  -  8 

2.5440D  -  8 

9.5736D  -  9 

14 

2.1389D  -  7 

9.7678D  -  8 

2.4503D  —  8 

1.3361D-8 

4.4863D  -  9 

15 

1.3910D  -  7 

6.0849D  -  8 

1.4005D  -  8 

7.3156D  -  9 

2.3041D  -  9 

HI 

3.9037D  +  0 

4.6895D  +  0 

6.8476D  +  0 

8.3136D  +  0 

1.2344D  +  1 
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TABLE  4.3.  The  Relative  Errors  in  Energy  Norm  for  a  =  m(3/2),m  =  2/3, 1, 

when  /  =  1. 


p  a 

1 

■H 

DOF 

L 

1 

0.6020654757148309D  4-  00 

0.5706890379221034D  +  00 

2 

0.1308957854381257D  +  00 

0.1009916680045339D  +  00 

3 

0.9925866863851720D  -  01 

0.8092169418184365D  -  01 

4 

0.42033661 19963193D  -01 

0.6942004795074728D  -  02 

100 

5 

0.3198798981482650D  -  01 

0.2059652008864891D  -  02 

152 

6 

0. 25638 98916834782D  -  01 

0.5410328220679518D  -  03 

216 

0.2124137148103968D  -  01 

0.1410620291464109D  -  03 

292 

0.1802105654621925D  -  01 

0.3647987063804444D  -  04 

380 

m 

■■1 

1 

0.6140576157783938D  +  00 

0.7227916745356786D  +  00 

KM 

2 

0. 1 153904981 172256D  +  00 

0.1685768303941887D  +  00 

pip* 

3 

0.8432233140432240D  -  01 

0.9410470781118301D  —  01 

pIpH 

4 

0.8046799392682128D  -  02 

0.28859981629781 18D  -01 

5 

0.2159258500150655D  -  02 

0.1095935137932326D  -  01 

152 

6 

0.5484874342364825D  -  03 

0.3570789379357526D  -  02 

216 

0. 1413509051434937D  -  03 

0.1052197244808133D  -  02 

292 

0.3612904449632896D  -  04 
...  _  „  -  .  . . 

0.30775796203461 13D  -  03 

380 

Finally  let  us  remark  that  if  /  ^  0  then  the  computation  of  the  elemental  load  vector  is 
influenced  by  the  auxiliary  mapping.  For  example,  if  the  conformal  mapping  ipa  :  £a  =  r 
was  used  as  in  example  3.1  then 

has  to  be  used  for  the  elemental  load  vectors  on  the  elements  in 

In  this  section  we  are  concerned  only  with  the  corner  singularities,  but  our  idea  can  be 
extended  to  deal  with  some  other  types  of  singularity  as  well.  For  example,  the  boundary 
singularity,  known  as  the  problem  of  Motz  in  the  literatures(see  [20],  [25],  [27]),  can  also 
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5.  The  Method  of  Auxiliary  Mapping  for  the  General  Elliptic  Differential 

Equation  of  Second  Order. 

In  the  previous  sections,  we  have  applied  the  auxiliary  mapping  technique  to  the  Laplace 
operator.  In  this  section  we  will  show  that  this  technique  can  also  be  used  for  general  elliptic 
boundary  value  problems.  In  the  case  of  Laplace  operator,  Corollary  2.1  guarantees  the 
simpleness  of  computing  the  elemental  stiffness  matrices.  However  in  the  case  of  general 
elliptic  operator  we  need  to  use  lemma  2.1  for  computing  elemental  stiffness  matrices  and 
hence  the  routines  for  differential  equations  with  nonconstant  coefficients  should  be  used 
for  the  computation  in  general.  Even  if  the  given  problem  has  constant  coefficients,  the 
coefficients  we  are  concerned  in  the  auxiliary  mapping  technique  are  polynomials  in  sin  x 
and  cos  x .  This  fact  has  to  be  taken  in  account  because  quadrature  formula  are  used  in  the 
computation  of  the  stiffness  matrices.  Nevertheless  they  lead  to  uniformly  elliptic  operator. 
Hence  for  /  =  0  one  can  achieve  once  again  the  exponential  rate  of  convergence  as  before. 
Similar  situation  also  occurs  for  /  ^  0. 

We  now  describe  the  mapping  technique  for  the  general  case  with  the  mesh  shown  in 
Fig.  3.1.  This  goes  as  follows: 

(Step  a)  Select  parameter  a  for  the  mapping 

(Step  b)  Map  elements  C7,...,ei2  onto  e*,...,ej2  by  the  mapping  (v?“)-1  and  construct  the 
coefficients  of  the  transformed  elliptic  operator  and  bilinear  form:  Since  (<^°)_1(r)  = 
zlla  we  have 


M  = 


cos(l  —  a)r/> 
sin(l  —  a)V> 


—  sin(l  —  a)V» 
cos(l  —  a)V> 


and  by  lemma  2.1  the  coefficient  qX]  of  the  transformed  elliptic  operator  can  be  com¬ 
puted  by  the  following  formula 


(5.1) 


1 


t  —  ( 1  —  q)0 

9n  =  du  cos2*  +  a22  sin2 t  -  (d2J  +  a12)sin*cos* 
qi2  =  (ai  i  —  a22 )  sin  t  cos  *  —  a2 1  sin2  *  +  dj 2  cos2  t 
q2i  =  (an  —  a22)sin  f  cos*  —  di2  sin2 1  +  a2i  cos2  t 
922  =  dn  sh*2  t  +  a22  cos2  t  +  (di2  -I-  a2i )  sin  t  cos  t 
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where  atJ  =  aXJ  o  ^>a  and  (p,  ip)  is  the  polar  coordiates  in  the  £-plane. 

(Step  c)  Taking  care  of  numerical  integration,  compute  the  elemental  stiffness  matrix  on  e* 
with  the  coefficients  qi;(£)  and  the  elemental  load  vector  on  e*  with  the  right  hand 

side  /  =  •  (/(/>,</>)) 

(Step  d)  Determine  the  finite  element  solution  ufe  by  solving  the  assembled  global  system. 

The  solution  at  a  point  x  £  ej  can  be  obtained  by  evaluating  the  solution  on  e*  at 


It  is  worth  to  note  the  following  observation  about  the  case  when  the  coefficient  atJ 
are  constants:  suppose  Aj  and  A2  are  the  eigenvalues  of  the  symmetric  matrix  A  =  [atJ] 
then  there  is  an  orthogonal  matrix  P  such  that  PAPT  —  diag  (Ai,A2).  If  'P  denotes  the 
mapping  defined  by 

*{z’y)={-kx’7r>y) 

then  *P  maps  the  elliptic  subdomain 

2  2 

JV  =  {(x,v) :  <  1}  nn 

Al  A2 


onto  the  circular  subdomain  =  {z  :  \z\  <  1}  fl  fi.  Let  us  denote  P  1{K)  and 
(¥,a)_1(^Ro)  by  K  and  Q,*Rq  respectively.  That  is, 


K  K  SIr0  (^‘  fit 


Ro- 


Then  by  an  argument  similar  to  that  of  lemma  2.1  we  have 


(5.2) 


Vu  •  Vvd( 
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In  particular,  if  ai2  =  021  =  0,  we  have 


(5-3) 


du  dv 
dxi  dxi 


+  a2  2 


du  du 
dx2  dx2 


)dx  = 


Vu  •  Vudf 


and  the  computation  of  elemental  stiffness  matrix  and  load  vector  on  e*  in  f 1*Rg  is  as  simple 
as  the  case  of  Laplace  operator.  We  can  also  use  here  the  triangular  elements  with  elliptic 
arcs  which  are  allowed  in  PROBE.  The  mapping  Technique  together  with  this  observation 
is  leading  to  another  approach  to  deal  with  interface  singularities.  Detailed  argument  of 
this  will  be  presented  elsewhere. 

In  similar  way,  we  can  treat  the  plane  elasticity  problem.  For  example,  in  the  case  of 
crack  where  the  sides  axe  not  loaded  the  solution  has  a  simple  form  and  the  auxiliary 
mapping  technique  with  the  conformal  mapping  2  =  will  lead  to  the  rate  N~a  of 
convergence  analogously  as  before.  Unfortunately  in  contrast  to  the  case  of  the  Laplace 
operator  the  local  stiffness  matrices  has  to  be  computed  via  nonconstant  coefficients  of  the 
form  similar  to  (5.1). 


6.  Elliptic  Problems  on  Unbounded  Domains. 

In  a  similar  manner  to  the  previous  sections  the  method  of  auxiliary  mappings  can  deal 
the  problems  on  unbounded  domains.  Let  us  address  this  with  a  model  problem  on  an 
unbounded  domain  ft  with  bounded  boundary  T  (this  is  ,  ft  is  the  outside  of  a  bounded 
domain  enclosed  by  a  simple  closed  curve  T). 


(6.1a) 

1 

l> 

e 

II 

in  ft, 

(6.1b) 

u  =  0 

on  VD 

du 

(6.1c) 

on  r  n 
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where  T d  and  Tjv  are  unions  of  analytic  arcs  as  section  3  and  T  is  a  simple  closed  curve. 
For  simplicity  we  will  assume  that  T £>  ^  0  and  the  orgin  is  inside  of  T. 

Let  be  the  conformal  mapping  £  =  -j,  and  let 


*)  =  {u  £  Le(il)  :{jj  u2|J($)|  +  //  Vu  •  Vu)  <  oo}, 
WlD{ ft)  =  {u  €  W\n,  tf)  :  u  =  0  on  TD) 


Suppose 


JL 


i 


nm  WW<00' 


then  there  exists  a  unigue  u  €  such  that 


I j  Vu  ■  Vvdx  =  If  fvdx  4-  /  gvds,  for  all  v  € 

J  J  n  -/r* 


Suppose  the  mesh  on  SI  is  as  shown  in  Fig.  6.1.  Then  by  the  conformal  mapping  'F 
the  infinite  elements  e°°  will  be  mapped  to  the  curvilinear  trangular  elements  e*.  Now 
we  can  proceed  as  before.  For  example,  suppose  /  =  0  and  on  T  =  dQ.  we  prescribed 
the  same  boundary  conditions  as  section  3  and  the  mesh  on  ft  is  as  shown  in  Fig.  6.1. 
Then  combining  the  auxiliary  mapping  technigue  in  the  neighborhood  of  every  vertex  and 
the  auxiliary  mapping  'P  for  the  infinite  elements  we  can  achieve  an  exponential  rate  of 
convergence  as  before,  i.e. 

||wo  —  Up||wi(n,'j')  <  Ce  7\/^ 

Let  us  show  one  numerical  example. 

Example  6.1.  Let  us  consider  Laplace’s  equation  on  ft  =  {z  :  |z|  >  1/2}.  Suppose 
£1Rq  =  {2  :  \z\  <  l}  ,  the  mesh  is  as  shown  in  Fig.  6.2.  and  we  impose  the  nonhomoge- 
neous  Dirichret  boundary  condition  on  T  by  the  true  solution  u(x,  y )  =  x/(x2  +  y 2).  Then 
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the  relative  error  in  energy  norm  is  as  shown  in  Table  6.1  and  we  cam  see  once  again  the 
exponential  rate  of  convergence  in  Fig.  6.3. 


Table  6.1.  The  Total  Strain  Energy  and  The  Relative  Error 
for  Laplace’s  Equation  on  the  Unbounded  Domain. 


p 

Total  Energy 

Relative  Error  in  Energy  Norm 

DOF 

1 

11.16236197243942 

0.3342565747335926D  +  00 

9 

2 

12.43592792490380 

0.1018837543895582D  + 00 

29 

3 

12.56842346856016 

0. 1278127328426342D  -  01 

53 

4 

12.56658836901951 

0.4162735313104689D  -  02 

89 

5 

12.56637844630835 

0.7894597600800805D  -  03 

137 

6 

12.56637082991567 

0. 130971 1457162764D  -  03 

189 

12.56637061959958 

0.2042102697231802D  -  04 

269 

12.56637061447803 

0.3075441953279617D  -  05 

353 

Similar  approach  can  also  be  applied  for  the  conformal  mapping  £  =  1  l(za)  and  also  for 
the  general  differential  equation  and  elasticity  problems.  It  is  essential  that  the  solution 
has  finite  energy  and  that  we  consider  the  behaviour  of  the  solution  at  oo.  So  far  we 
mentioned  only  the  case  when  the  boundary  is  bounded  only.  By  the  same  approach  we 
can  deal  with  infinite  domains  such  as  half  plane  etc.  We  memtion  that  the  implemention 
of  the  infinite  elements  is  exactly  ia  simple  as  for  ther  finite  elements  we  disscused  in 
previous  sections.  Many  times  we  are  dealing  with  the  case  that  the  solution  does  not 
belong  to  the  space  Hl(Cl)  because  of  the  behavour  at  oo.  Typical  case  is  when  on  the 
bounadry  F  we  prescribed  the  following  singular  boundary  condition: 


du 

dn 


=  9, 


and  /  =  0. 
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Then  the  solution  at  oc  is  of  form  u  «  C  log  r  and  we  have  to  enrich  the  space  of  elemen¬ 
tal  shape  functions  by  this  function  analoguously  as  we  disscussed  earlier.  Of  course  we 
have  to  adjust  the  approach  to  the  fact  that  we  deal  with  one  shape  function  with  infinite 
energy.  This  can  be  dealt  in  different  manner. 


7.  Conclusion 

The  h-p  version  allows  a  very  effective  and  simple  treatment  of  the  singularities  of 
the  solution  caused  by  the  comers  of  domains  and  by  the  unboundedness  of  domains  in 
R2 .  This  approach  can  be  combined  with  elemental  enrichment  procedure  and  the  mesh 
refinement  procedure. 

By  using  this  approach  one  can  obtain  a  higher  rate  of  convergence  in  almost  all  appli¬ 
cations.  Prom  the  view-point  of  implementation,  the  easiest  and  the  cheapest  approach 
is  the  auxiliary  mapping  technique  with  possibly  nonunifom  polynomial  degrees  over 
the  elements. 
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